#coding=utf-8
import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate

def g(epsilon):
    if epsilon < 0.001:
        return epsilon**2 * np.exp(epsilon)
    else:
        return epsilon**4 * np.exp(-epsilon) / (1 - np.exp(-epsilon) )**2;
'''
x = np.arange(0.1,10,0.1)
y = [g(t) for t in x]
plt.plot(x,y)
plt.show()
'''
def Cv(T):
    #print("1.0/T = ", 1.0/T)
    y = 9 * 8.31 * T**3 * integrate.quad( g, 0, 1.0/T )[0]
    return y

T = np.arange(1e-5, 1, 0.01 )
C = [Cv(t) for t in T]
plt.plot(T,C)
plt.show()